function [P_l_up,T_l_up,rou_l_up,Ma_l_up,a_l_up,P_l_down,T_l_down,rou_l_down,Ma_l_down,a_l_down]...
    =trim(kk,P_air,T_air,rou_air,Ma,coef,aero_up_pre,aero_down_pre,stru_up_rigid,stru_down_rigid,stru_grid,aero_up_con,aero_down_con,aerobox,Npath)

for i=1:size(stru_up_rigid,1)
    stru_up_rigid(i,2:4)=stru_grid(stru_grid(:,1)==stru_up_rigid(i,1),2:4);
end
for i=1:size(stru_down_rigid,1)
    stru_down_rigid(i,2:4)=stru_grid(stru_grid(:,1)==stru_down_rigid(i,1),2:4);
end

[~,WAD_up_pre_rigid]=interpolate(stru_up_rigid(:,2),stru_up_rigid(:,3),stru_up_rigid(:,4),aero_up_pre(:,2),aero_up_pre(:,3));
[~,WAD_down_pre_rigid]=interpolate(stru_down_rigid(:,2),stru_down_rigid(:,3),stru_down_rigid(:,4),aero_down_pre(:,2),aero_down_pre(:,3));

try
    load temp\A.mat A;
catch
    A=[0;0];
end

l=1;
alpha=A(1);
rudder=A(2);
e=0.000005;
dalpha_temp=0.05; drudder_temp=0.05; f_temp=10^8; m_temp=10^8; df_temp=10^10; dm_temp=10^10; f=[0,0,0]; m=[0,0,0];

if isempty(gcp('nocreate'))
    parpool(3);
end
dir=["101o","101f","101m"];
for ii=dir
    i=char(ii);
    mkdir(i);
    eval(['copyfile 101.bdf ',i,'\',i,'.bdf']);
    eval(['copyfile grid_101.bdf ',i,'\grid_101.bdf']);
    eval(['copyfile bar.bdf ',i,'\bar.bdf']);
    eval(['copyfile tria.bdf ',i,'\tria.bdf']);
    eval(['copyfile quad.bdf ',i,'\quad.bdf']);
    eval(['copyfile penta.bdf ',i,'\penta.bdf']);
    eval(['copyfile hexa.bdf ',i,'\hexa.bdf']);
    eval(['copyfile spc.bdf ',i,'\spc.bdf']);
    eval(['copyfile rbe.bdf ',i,'\rbe.bdf']);
    eval(['copyfile rbe2.bdf ',i,'\rbe2.bdf']);
    eval(['copyfile property.bdf ',i,'\property.bdf']);
    eval(['copyfile material_101.bdf ',i,'\material_101.bdf']);
    try
        eval(['copyfile 103.pch ',i,'\103.pch']);
    catch
    end
end

while true
    input_alpha=[alpha,alpha+e,alpha]; input_rudder=[rudder,rudder,rudder+e];
    for i=1:3
        [P_l_up,~,~,~,~,P_l_down,~,~,~,~]...
            =shock_expansion(kk,P_air,T_air,rou_air,Ma,input_alpha(i),input_rudder(i),aero_up_pre,aero_down_pre,WAD_up_pre_rigid,WAD_down_pre_rigid);
        P_l_up=1/coef*P_l_up; P_l_down=1/coef*P_l_down;
        create_101_P(P_l_up,P_l_down,aero_up_con,aero_down_con,aerobox);
        dir_path=char(dir(i));
        eval(['copyfile pressure.bdf ',dir_path,'\pressure.bdf']);
        eval(['copyfile load.bdf ',dir_path,'\load.bdf']);
    end
    delete('pressure.bdf');
    parfor i=1:3
        dir_path=char(dir(i));
        run_101(Npath,dir_path);
        delete(['101',dir_path(end),'.pch']);
        [f(i),m(i)]=read_101_force(dir_path);
        delete(['101',dir_path(end),'.f06']);
    end

    if f(1)==f_temp&&m(1)==m_temp
        dalpha_temp=0.001;
        drudder_temp=0.001;
        l=0.5;
    end
    df=f(1)-f_temp; dm=m(1)-m_temp; f_temp=f(1); m_temp=m(1);
    if abs(df)>0.9*abs(df_temp)||abs(dm)>0.9*abs(dm_temp)
        l=0.8*l;
    end
    df_temp=df; dm_temp=dm;

    disp(f(1)/(50*10^3))
    disp(m(1)/(50*10^6))
    if abs(f(1))<50*10^3&&abs(m(1))<50*10^6
        break
    end

    dalpha_temp=0.8*dalpha_temp;
    drudder_temp=0.8*drudder_temp;
    d1=(f(2)-f(1))/e; d2=(f(3)-f(1))/e; d3=(m(2)-m(1))/e; d4=(m(3)-m(1))/e;
    if d2==0&&d4==0
        A(1)=alpha-l*f(1)/d1;
        A(2)=rudder+0.000001;
    elseif (d1/d2)==(d3/d4)
        A(2)=rudder-l*f(1)/d2;
    else
        D=[d1 d2;d3 d4]\[f(1);m(1)];
        if dalpha_temp<abs(D(1))
            D(1)=0.6*dalpha_temp*D(1)/abs(D(1));
        end
        if drudder_temp<abs(D(2))
            D(2)=0.6*dalpha_temp*D(2)/abs(D(2));
        end
        A(1)=alpha-l*D(1);
        A(2)=rudder-l*D(2);
    end
    A(1)=min(A(1),0.09); A(1)=max(A(1),-0.09);
    A(2)=min(A(2),0.17); A(2)=max(A(2),-0.17);
    alpha=A(1);
    rudder=A(2);
    save temp\A.mat A;
end
delete(gcp('nocreate'));
rmdir('101o','s'); rmdir('101f','s'); rmdir('101m','s');
[P_l_up,T_l_up,rou_l_up,Ma_l_up,a_l_up,P_l_down,T_l_down,rou_l_down,Ma_l_down,a_l_down]...
    =shock_expansion(kk,P_air,T_air,rou_air,Ma,alpha,rudder,aero_up_pre,aero_down_pre,WAD_up_pre_rigid,WAD_down_pre_rigid);
end